The Challenge & Limits of SDMs

• Critical Challenge: Predicting species’ responses to global environmental change is vital for conservation and management.

• Usual Tool: Species Distribution Models (SDMs) have been the workhorse, correlating occurrences with environmental variables.

• Limitations of SDMs:

  • Struggle to predict responses under novel future conditions (Pagel and Schurr 2012).
  • Lack Mechanism: They do not explicitly model the underlying biological processes.
  • Equilibrium Assumption: Often violated (Guisan and Thuiller 2005).

Dynamic Range Models (DRMs): A Mechanistic Approach

• DRMs explicitly incorporate demographic processes that drive range dynamics (Pagel and Schurr 2012).

• Allows linking environmental drivers directly to specific processes.

• Potential for more robust forecasting under novel conditions.

• Why are they rarely in practice? Despite conceptual appeal, DRMs have been underutilized due to their complexity and computational challenges to fit these models.

Introducing drmr: Making DRMs Accessible

• Goal: Bridge the gap between DRM potential and practical application.

• Key Features:

  • Makes age-structured DRMs readily available in a Bayesian framework.
  • Leverages Stan via cmdstanr for efficient fitting (Gabry et al. 2024).
  • User-friendly interface.
  • Easily relate environmental drivers to specific recruitment and survival.

Setup & Notation

• \(Y_{a, t, i}\): Unobserved density of individuals of age \(a\), time \(t\) and site \(i\).

• \(Y_{t, i} = \sum_{a} Y_{a, t, i}\): Observed density of individuals (all ages) at time \(t\) and site \(i\).

• \(\lambda_{a, t, i}\): Expected age-specific density.

  • Biological processes (recruitment, survival, and movement) are encoded through these parameters.

• \(\mu_{t, i} = \sum_{a} \lambda_{a, t, i}\): Expected total density.

• \(\rho_{t, i} = \mathrm{P}(Y_{t, i} = 0)\): Probability of absence.

Modeling assumption: Zero-augmented probability density function (pdf)

\[ f(y_{t, i} \mid \mu_{t, i}, \phi, \rho_{t, i}) = \begin{cases} \rho_{t, i}, & \text{ if } y_{t, i} = 0, \\ (1 - \rho_{t, i}) g(y_{t, i} \mid \mu_{t, i}, \phi), & \text{ if } y_{t, i} > 0. \end{cases} \]

• \(y_{t, i}\) is a realization of \(Y_{t, i}\).

• \(g(\cdot \mid \mu_{t, i}, \phi)\) is the pdf of a continuous probability distribution.

• \(\phi\) is a nuisance parameter.

Let’s look at some model assumptions!

Recruitment

graph TD
    xr["$$\Large  \mathbf{x}^{(r)}_{t, i}$$"] --> lambda1(("$$\Large \lambda_{1, t, i}$$"));
    betar(("$$\Large \boldsymbol{\beta}_r$$")) --> lambda1;
    zr(("$$\Large z^{(r)}_{t, i}$$")) --> lambda1;

Survival

graph TD
    subgraph "Time t-1"
        xs["$$\Large \mathbf{x}^{(s)}_{t - 1, i}$$"] --> s_prev;
        betas(("$$\Large  \boldsymbol{\beta}_s$$")) --> s_prev;
        zs(("$$\Large  z^{(s)}_{t - 1, i}$$")) --> s_prev;
        f["$$\Large  f_{a - 1, t - 1}$$"] --> s_prev;
        lambda_prev(("$$\Large \lambda_{a - 1, t - 1, i}$$"))
        s_prev(("$$\Large s_{a - 1, t - 1, i}$$"))
    end
        
    lambda_prev --> lambda_a(("$$\Large \lambda_{a, t, i}$$"));
    s_prev --> lambda_a;

Process and Observation error

graph TD
    %% Mean calculation
    lambda_a(("$$\Large \lambda_{a, t, i}$$"));
    lambda_a -- $$\Large \sum_a \lambda_{a, t, i}$$ --> mu(("$$\Large \mu_{t, i}$$"));

subgraph "P(Y = 0)"
    xt["$$\Large \mathbf{x}^{(p)}_{t, i}$$"] --> rho(("$$\Large \rho_{t, i}$$"));
    betat(("$$\Large \boldsymbol{\beta}_p$$")) --> rho;
end

%% Final Observation
    mu --> Y{"$$\Large Y_{t, i}$$"};
    %% Other components
    phi(("$$\Large \phi$$")) --> Y;
    rho --> Y;

Demographic processes

graph TD
    subgraph "Demographic processes"
        %% Inputs and Parameters for lambda1
        xr["$$\mathbf{x}^{(r)}_{t, i}$$"] --> lambda1(("$$\lambda_{1, t, i}$$"));
        betar(("$$\boldsymbol{\beta}_r$$")) --> lambda1;
        zr(("$$z^{(r)}_{t, i}$$")) --> lambda1;

        %% Recursive part for lambda_a
        subgraph "Time t-1"
            xs["$$\mathbf{x}^{(s)}_{t - 1, i}$$"] --> s_prev;
            betas(("$$\boldsymbol{\beta}_s$$")) --> s_prev;
            zs(("$$z^{(s)}_{t - 1, i}$$")) --> s_prev;
            f["$$f_{a - 1, t - 1}$$"] --> s_prev;
            lambda_prev(("$$\lambda_{a - 1, t - 1, i}$$"))
            s_prev(("$$s_{a - 1, t - 1, i}$$"))
        end
        
        lambda_prev --> lambda_a(("$$\lambda_{a, t, i}$$"));
        s_prev --> lambda_a;
    end

    %% Mean calculation
    lambda_a --> mu(("$$\mu_{t, i}$$"));

subgraph "Zero-augmentation"
        xt["$$\mathbf{x}^{(p)}_{t, i}$$"] --> rho(("$$\rho_{t, i}$$"));
        betat(("$$\boldsymbol{\beta}_p$$")) --> rho;
    end

    %% Final Observation
    mu --> Y{"$$Y_{t, i}$$"};

    %% Other components
    phi(("$$\phi$$")) --> Y;
    rho --> Y;

Code & workflow

Summer Flounder Dataset

• An example analysis uses Summer flounder (Paralichthys dentatus) data from 1982-2016 NOAA bottom trawl surveys (Fredston et al. 2025) to illustrate the package’s features.

• The data spans the US Atlantic coast (Cape Hatteras, NC to the Canada/Maine border) and was aggregated from individual hauls into 10 latitudinal patches with varying areas.

• Response variable: Density (count per unit area).

• Environmental drivers: SST and SBT.

Models fitted to the data

Code for DRM (rec)

drm_rec <-
  fit_drm(.data = dat_train,
          y_col = "dens", ## response variable: density
          time_col = "year", ## vector of time points
          site_col = "patch",
          family = "gamma", ## other options: lognormal, loglogistic
          seed = 202505,
          formula_zero = ~ 1 + c_hauls,
          formula_rec = ~ 1 + c_sst + I(c_sst * c_sst),
          formula_surv = ~ 1,
          f_mort = f_train,
          n_ages = NROW(f_train),
          adj_mat = adj_mat, ## A matrix for movement routine
          ages_movement = c(0, 0,
                            rep(1, 12),
                            0, 0), ## ages allowed to move
          .toggles = list(ar_re = "rec", ## other options: "surv", "dens"
                          movement = 1,
                          est_surv = 1,
                          est_init = 0,
                          minit = 1),
          .priors = list(pr_phi_a = 1, pr_phi_b = .1,
                         pr_alpha_a = 4.2, pr_alpha_b = 5.8,
                         pr_zeta_a = 7, pr_zeta_b = 3))

Code for DRM (surv)

drm_surv <-
  fit_drm(.data = dat_train,
          y_col = "dens", ## response variable: density
          time_col = "year", ## vector of time points
          site_col = "patch",
          family = "gamma", ## other options: lognormal, loglogistic
          seed = 202505,
          formula_zero = ~ 1 + c_hauls,
          formula_rec = ~ 1,
          formula_surv = ~ 1 + c_sbt + I(c_sst * c_sbt),
          f_mort = f_train,
          n_ages = NROW(f_train),
          adj_mat = adj_mat, ## A matrix for movement routine
          ages_movement = c(0, 0,
                            rep(1, 12),
                            0, 0), ## ages allowed to move
          .toggles = list(ar_re = "rec", ## other options: "surv", "dens"
                          movement = 1,
                          est_surv = 1,
                          est_init = 0,
                          minit = 1),
          .priors = list(pr_phi_a = 1, pr_phi_b = .1,
                         pr_alpha_a = 4.2, pr_alpha_b = 5.8,
                         pr_zeta_a = 7, pr_zeta_b = 3))

Effect of environmental variables

Forecasting

• The predict_drm function requires at least three arguments:
  • drm: Output of a fit_drm call
  • new_data: A data.frame where we seek to obtain predictions
  • seed: For reproducibility

forecast_rec <- predict_drm(drm = drm_rec,                         # Required
                            new_data = dat_test,                   # Required
                            seed = 152,                            # Required
                            f_test = f_test,                       # Optional
                            past_data = filter(dat_train,          # Optional
                                               year == max(year)))

• When \(f_{a, t}\) is used at the time of model fit, we also need to provide f_test to predict_drm.

• past_data is necessary when survival is a function of environmental variables.

Forecasts visualization & comparison

Highlights

• The drmr substantially lowers the barrier for ecologists to use the DRM in their applications.

• The code is easy to use and takes advantage of what has been developed for Stan: visualization, diagnostic tools, and estimation.

• drmr allows for empirically testing which processes are more important to predict the distribution of a species.

• The more complex a model is, the more (and better) data we need to be able to estimate those relationships.

Future work

• Include population dynamic models that explicitly relate adult density to recruitment (e.g., Ricker, Belverton-Holt)

• GAM-like non-parametric relationships between processes and the environment

• Support for length-composition data

• More realistic movement routines

Acknowledgments